global SSDIMed "/disk/agedisk4/medicare.work/miller-DUA50377/proj_ssdi"
* Settings
version 16
do "$SSDIMed/scripts/_auxiliary/_project_settings.do"
do "$SSDIMed/scripts/_auxiliary/_project_grstyle.do"

* All required Stata modules are available in the /_packages and /_auxiliary folders
cap adopath - "/home/site/etc/stata/ado.nber"
cap adopath - "/home/site/etc/stata/ado.nber/updates"
adopath ++ "$SSDIMed/scripts/_packages"
adopath ++ "$SSDIMed/scripts/_auxiliary"
grstyle init 



*0) BOOTSTRAP SAMPLE
/*
set seed 20170526
cap mkdir "$SSDIMed/data/analysis/bootstrap"
local bsnum 500

forval n=1/`bsnum'{
use fipscounty_firstnm  covstart_month using "$SSDIMed/data/analysis/county-month-age_entry_sample-main.dta", clear 
rename (fipscounty_firstnm ) (fipscounty_init)
gduplicates drop _all, force 
bsample 
rename * *`n'
gen clustervar`n'=_n
gen counter=_n 
save "$SSDIMed/data/analysis/bootstrap/sample`n'.dta", replace
}

use "$SSDIMed/data/analysis/bootstrap/sample1.dta", clear 
forval n=2/`bsnum'{
qui merge 1:1 counter using "$SSDIMed/data/analysis/bootstrap/sample`n'", nogen 
}
save "$SSDIMed/data/analysis/bootstrap/samples.dta", replace 
*/

*************************************************************************************************
*1) BOOTSTRAP FOR ENTRY
/*
foreach s in 0 1{
local sample 5152
use "$SSDIMed/data/analysis/county-month-agegender_entry_sample-main.dta", clear  
keep if male_init==`s'
keep if inrange(age_year_covstart_fill, 51, 52)
gen age52 = (age_year_covstart_fill == 52)
sum unemp_rate_county_atapp [aw = pop_19_61_atapp]
gen UR=unemp_rate_county_atapp-`r(mean)'
gen age52xUR=age52*UR
***verify same as current results
local controls_02 i.fipscounty_firstnm_g
local y incidence_pop_agegender_atapp
local w pop_agegender_atapp
local clust county_mofd
local ctrl "02"
reghdfe `y' age52 UR age52xUR [aw = `w'], abs(`controls_`ctrl'') cluster(`clust')

cd "$SSDIMed/results/estimates/x-age52xUR/"
cap postclose entry_unempatappx52
postfile entry_unempatappx52 intercept UR age52  age52xUR sample using "bootstrap_incidence_male`s'.dta", replace 
local bsnum 500
forval n=1/`bsnum'{
di "`n'"
qui{
use "$SSDIMed/data/analysis/county-month-agegender_entry_sample-main.dta", clear  
keep if male_init==`s'
keep if inrange(age_year_covstart_fill, 51, 52)
gen age52 = (age_year_covstart_fill == 52)
sum unemp_rate_county_atapp [aw = pop_19_61_atapp]
gen UR=unemp_rate_county_atapp-`r(mean)'
gen age52xUR=age52*UR
rename (fipscounty_firstnm covstart_month ) (fipscounty_init`n' covstart_month`n')
joinby fipscounty_init`n' covstart_month`n' using "$SSDIMed/data/analysis/bootstrap/samples.dta"
rename (fipscounty_init`n' covstart_month`n') (fipscounty_init covstart_month ) 
local y incidence_pop_agegender_atapp 
local clust county_mofd
local w pop_agegender_atapp
reghdfe `y' age52 UR age52xUR [w=`w'], abs(fipscounty_init)  cluster(`clust')
local intercept=_b[_cons]
local UR=_b[UR]
local age52=_b[age52]
local age52xUR=_b[age52xUR]
post entry_unempatappx52 (`intercept') (`UR') (`age52') (`age52xUR') (`n') 
}
}
postclose entry_unempatappx52
}
*/



*************************************************************************************************
*2) BOOTSTRAP FOR SPENDING 
/*
foreach s in 0 {
cd "$SSDIMed/results/estimates/x-age52xUR"
postfile spending_unempatappX52 intercept age52 UR age52xUR sample using "bootstrap_spending_male`s'.dta", replace 
forval n=1/500{
di "`n'"
qui{
use "$SSDIMed/data/analysis/bene-year_outcomes_sample-main.dta",  clear
keep if inrange(age_year_covstart_fill, 51, 52)
keep if male_init==`s'
rename (fipscounty_firstnm covstart_month ) (fipscounty_init`n' covstart_month`n')
joinby fipscounty_init`n' covstart_month`n' using "$SSDIMed/data/analysis/bootstrap/samples.dta"
gegen fipscounty_init=group(fipscounty_init`n')
gen age52 = (age_year_covstart_fill == 52)
sum unemp_rate_county_atapp 
gen UR=unemp_rate_county_atapp-`r(mean)'
gen age52xUR=age52*UR
local clust county_mofd
local y tot_pmt 
local controls_02 i.years_since_covstart##i.fipscounty_init
reghdfe `y' age52 UR age52xUR, abs(`controls_02') 
local intercept=_b[_cons]
local UR=_b[UR]
local age52=_b[age52]
local age52xUR=_b[age52xUR]
post spending_unempatappX52 (`intercept') (`age52') (`UR') (`age52xUR') (`n') 
}
}
postclose spending_unempatappX52
}
*/

*************************************************************************************************
*3) bring together true and bootstrap results for alphas and betas  
/*
foreach s in 0 1{
clear 
use "$SSDIMed/results/estimates/robustness/Column1robustnessB_male`s'.dta", 
drop if var=="" | var=="controls"
replace var="intercept" if var=="_cons" 
keep var coef1 
gen i=1
reshape wide coef1, i(i) j(var) string 
rename coef1* *
drop i 
gen sample=0 
append using "$SSDIMed/results/estimates/x-age52xUR/bootstrap_incidence_male`s'.dta",
compress _all 
rename (intercept age52 UR age52xUR ) (alpha alpha52 alphaUR alpha52XUR)
save "$SSDIMed/results/estimates/modelparameterization/alphas_fipscounty_init_male`s'.dta", replace 
list in 1/10
sum 



clear 
use "$SSDIMed/results/estimates/robustness/Column2robustnessB_male`s'.dta", 
drop if var=="" | var=="controls"
drop in 8
drop in 6
drop in 4
drop in 2
replace var="intercept" if var=="_cons" 
keep var coef1  
gen i=1
reshape wide coef1, i(i) j(var) string 
rename coef1* *
drop i 
gen sample=0 
append using "$SSDIMed/results/estimates/x-age52xUR/bootstrap_spending_male`s'.dta",
compress _all 
rename (intercept age52 UR age52xUR ) (beta beta52 betaUR beta52XUR)
save "$SSDIMed/results/estimates/modelparameterization/betas_02_male`s'.dta", replace 
}
*/



*************************************************************************************************
*4) do calculations 
foreach s in 0 1{
use "$SSDIMed/results/estimates/modelparameterization/alphas_fipscounty_init_male`s'.dta", clear 
merge 1:1 sample using "$SSDIMed/results/estimates/modelparameterization/betas_02_male`s'.dta", nogen   


***benefit parameters
gen double m=2*beta52/alpha52 
gen double n=beta-(alpha*beta52/alpha52)
gen double mUR=2*(beta52+beta52XUR)/(alpha52+alpha52XUR) 
gen double nUR=beta+betaUR- (alpha+alphaUR)*(beta52+beta52XUR)/(alpha52+alpha52XUR) 
gen double mURminusm=2*(beta52+beta52XUR)/(alpha52+alpha52XUR) - 2*beta52/alpha52
gen double nURminusn=betaUR- (alpha+alphaUR)*(beta52+beta52XUR)/(alpha52+alpha52XUR) + (alpha*beta52/alpha52)
***elements of the figure
gen double x51=alpha 
gen double x51UR=alpha+alphaUR 
gen double x52=alpha+alpha52 
gen double x52UR=alpha+alpha52+alphaUR+alpha52XUR
gen double B_x51=n+m*x51 
gen double B_x51UR=nUR+mUR*x51UR 
gen double B_x52=n+m*x52
gen double B_x52UR=nUR+mUR*x52UR
***difference at x51UR and x52UR
gen double BU_x51=nUR+mUR*x51
gen double BU_x52=nUR+mUR*x52
***set up CI 
gen double B_200=n+m*200
gen double BUR_200=nUR+mUR*200 
gen double B_800=n+m*800
gen double BUR_800=nUR+mUR*800 
***cost parameters
*****COST PARAMETERS 
foreach DeltaC in -1000 -5000{
local posDeltaC=-1*`DeltaC'
gen double mC_51`posDeltaC'=(-`DeltaC' -m*alpha -n +mUR *(alpha+alphaUR) +nUR)/alphaUR
gen double nC_51`posDeltaC'=m*alpha +n- mC_51`posDeltaC'*alpha
gen double alphaCF_51_C`posDeltaC' =(`DeltaC' +nC_51`posDeltaC'-n)/(m-mC_51`posDeltaC')
gen double mC_52`posDeltaC'=(-`DeltaC'  +mUR*(alpha+alphaUR+alpha52+alpha52XUR) -m*(alpha+alpha52) +nUR - n)/(alphaUR+alpha52XUR)
gen double nC_52`posDeltaC'=m*(alpha+alpha52)+n - mC_52`posDeltaC'*(alpha+alpha52)
gen double alphaCF_52_C`posDeltaC' =(`DeltaC' +nC_52`posDeltaC'-n)/(m-mC_52`posDeltaC')
gen double percentCF_health`posDeltaC'=(x51UR -alphaCF_51_C`posDeltaC' + x52UR-alphaCF_52_C`posDeltaC')/(alphaUR+alphaUR+alpha52XUR)
}

foreach n in 51 52{
sum B_x`n' in 1 
local B_x`n'=`r(mean)'
gen shareBSBU_x`n'belowB_x`n'=(BU_x`n'<`B_x`n'')
replace shareBSBU_x`n'belowB_x`n'=. in 1 
sum BU_x`n' in 1 
local BU_x`n'=`r(mean)'
gen shareBSB_x`n'aboveBU_x`n'=(B_x`n'>`BU_x`n'')
replace shareBSB_x`n'aboveBU_x`n'=. in 1
}

*5) collapse and print table 
cap mkdir "$SSDIMed/results/temp"
cap mkdir "$SSDIMed/results/tables"
preserve 
keep in 1 
xpose, clear varname 
save "$SSDIMed/results/temp/tablecolumn`column'_value_male`s'.dta", replace 
restore 
preserve 
drop in 1 
gcollapse (mean) shareBS*, fast 
xpose, clear varname 
save "$SSDIMed/results/temp/tablecolumn`column'_pval_male`s'.dta", replace
restore 
preserve 
drop in 1 
gcollapse (sd) *, fast 
xpose, clear varname 
save "$SSDIMed/results/temp/tablecolumn`column'_sd_male`s'.dta", replace 
restore 
preserve 
drop in 1 
gcollapse (p5) *, fast 
xpose, clear varname 
save "$SSDIMed/results/temp/tablecolumn`column'_lb_male`s'.dta", replace 
restore 
preserve 
drop in 1 
gcollapse (p95) *, fast 
xpose, clear varname 
save "$SSDIMed/results/temp/tablecolumn`column'_ub_male`s'.dta", replace
clear 
gen stat=""
foreach var in value sd lb ub pval{
append using "$SSDIMed/results/temp/tablecolumn`column'_`var'_male`s'.dta"
replace stat="`var'" if stat==""
}
reshape wide v1, i(_varname) j(stat) string 
rename v1* * 
rename _varname var 
drop if var=="sample"
save "$SSDIMed/results/estimates/modelparameterization/tablecolumn`column'_male`s'.dta", replace 
restore 
}
*/



****OUTSHEET TO EXCEL 
foreach s in 0 1{
clear 
append using "$SSDIMed/results/estimates/modelparameterization/tablecolumn`column'_male`s'.dta",

***drop figure elements 
gen column="1" 
gen dropme=inlist(var, "x51", "x52",  "", "B_x51UR", "B_x52", "B_x52UR")
replace dropme=1 if inlist(var, "B_200", "BUR_200", "B_800", "BUR_800")
replace dropme=1 if inlist(var, "alpha", "alpha52", "alpha52XUR", "alphaUR")
replace dropme=1 if inlist(var, "beta", "beta52", "beta52XUR", "betaUR")
***keep only deltaC=-5000, except for Panel C stuff 
local posDeltaC 5000
replace dropme=1 if regexm(var, "500")
replace dropme=0 if regexm(var, "5000")
replace dropme=1 if regexm(var, "50000")
replace dropme=0 if regexm(var, "alphaCF")
replace dropme=1 if regexm(var, "hi") | regexm(var, "lo")
drop if dropme==1
drop ub lb dropme 
reshape wide value sd, i(var) j(column) string 
gen counter=1 if var=="m"
replace counter=2 if var=="n"
replace counter=3 if var=="mUR"
replace counter=4 if var=="nUR"
replace counter=5 if var=="mURminusm"
replace counter=6 if var=="nURminusn"
replace counter=6.25 if var=="BU_x51"
replace counter=6.5 if var=="B_x51"
replace counter=6.75 if var=="diff_x51"
replace counter=7 if var=="mC_51`posDeltaC'"
replace counter=8 if var=="nC_51`posDeltaC'"
replace counter=9 if var=="mC_52`posDeltaC'"
replace counter=10 if var=="nC_52`posDeltaC'"
replace counter=11 if var=="x51UR"
replace counter=12 if var=="alphaCF_51_C500"
replace counter=13 if var=="alphaCF_51_C5000"
replace counter=14 if var=="alphaCF_51_C50000"
replace counter=15 if var=="x52UR"
replace counter=16 if var=="alphaCF_52_C500"
replace counter=17 if var=="alphaCF_52_C5000"
replace counter=18 if var=="alphaCF_52_C50000"
sort counter 
drop counter 
order var value1 sd1 
export excel $SSDIMed/results/tables/modelparameterization_male`s'.xls, firstrow(variables) replace 
}
*/

foreach s in 0 1 {
  local DeltaC=-5000
  local posDeltaC=-1*`DeltaC'
  clear 
  append using "$SSDIMed/results/estimates/modelparameterization/tablecolumn`column'_male`s'.dta",
  drop sd 
  drop if (regexm(var, "500") & !regexm(var, "5000")) | regexm(var, "50000")
  keep if regexm(var, "x51") | regexm(var, "x51UR") | regexm(var, "x52") | regexm(var, "x52UR") | regexm(var, "00")
  drop if regexm(var, "C_") | regexm(var, "diff_")
  replace var=subinstr(var, "`posDeltaC'", "", 1) 
  gen x=.
  gen xlb=.
  gen xub=.
  foreach var in x51 x51UR x52 x52UR {
    replace x=value if var=="`var'"
    replace xlb=lb if var=="`var'"
    replace xub=ub if var=="`var'"
  }
  
  replace x=200 if var=="B_200" | var=="BUR_200"
  replace x=800 if var=="B_800" | var=="BUR_800"
  gen B=. 
  gen Blb=.
  gen Bub=.
  foreach var in B_x51 B_x51UR B_x52 B_x52UR B_200 BUR_200 B_800 BUR_800 {
    replace B=value if var=="`var'"
    replace Blb=lb if var=="`var'"
    replace Bub=ub if var=="`var'"
  }
  gen UR=regexm(var, "UR")
  replace var=subinstr(var, "B_", "", 1)
  replace var=subinstr(var, "BUR_", "", 1)
  keep if inlist(var, "x51", "x51UR", "x52", "x52UR", "200", "800")
  foreach var in B Blb Bub {
    replace `var'=`var'/1000
  }

    ***puts everything in one line 
  gcollapse x xlb xub B Blb Bub  , by(var UR ) 
  foreach var in x51 x51UR x52 x52UR {
    sum x if var=="`var'" 
    local `var': di %3.0f = `r(mean)'
  }
  
  * Override default graph sizes
  grstyle set size 8pt: text_option
  grstyle set size 9pt: axis_title
  grstyle set size 8pt: tick_label
  grstyle set size 2pt: tickgap
  
  * Arrows and lines
  local p 5
  grstyle set color gs8: p`p'
  grstyle set lpattern solid: p`p'
  grstyle set linewidth vvthin: p`p'
  
  * How many major ticks?
  * How much should range extend beyond ticks (as fraction of tick step size)
  local y_ticks 8
  local y_margin_lo 0.4
  local y_margin_hi 0.4
  
  * y1 tick lo, step size, format
  local y1_tick_lo 9
  local y1_step 1
  local y1_fmt "%9.1f"
  local y1_sign "$"
  local y1_unit " thousand"
  
  * Construct y labels and scales from parameters above
  forvalues y = 1/1 {
    local y`y'_tick_hi  = `y`y'_tick_lo' + `y`y'_step' * (`y_ticks' - 1)
    local y`y'_range_lo = `y`y'_tick_lo' - `y`y'_step' * `y_margin_lo'
    local y`y'_range_hi = `y`y'_tick_hi' + `y`y'_step' * `y_margin_hi'
    di "y`y'_tick_hi:  `y`y'_tick_hi'"
    di "y`y'_range_lo: `y`y'_range_lo'"
    di "y`y'_range_hi: `y`y'_range_hi'"
    
    local y`y'_label ylabel(`y`y'_tick_lo'(`y`y'_step')`y`y'_tick_hi', axis(`y') format("`y`y'_fmt'") noticks nolabels labgap(1pt)) 
    local y`y'_scale yscale(`yalt' range(`y`y'_range_lo' `y`y'_range_hi') axis(`y') noline) 
    di `"y`y'_label:  `y`y'_label'"'
    di `"y`y'_scale:  `y`y'_scale'"'
  }
  
  * x-axis settings
  local xlo 120
  local xhi 830
  local xscale xscale(range(`xlo' `xhi'))
  local x_axis xlabel(200(200)800) `xscale' xtitle("Entrants per million", margin(t=-.5))
  
  * style and color for y-axis title/labels
  local y1_title_style it
  local y1_title_color gs8
  local y1_label_style it
  local y1_label_color gs8
  local y1_label_size 8pt
  
  local text_style sf
  local text_color black
  
  * y-axis labels placed above grid rule (build manually)
  local y1_x `xlo'
  local y1_place ne
  local y1_just left
  local y1_width 25
  
  local y 1
  local y`y'_la
  foreach tick of numlist  `y`y'_tick_lo'(`y`y'_step')`y`y'_tick_hi' {
    if `tick' == `y`y'_tick_hi' {
      * white textbox to go behind label
      local y`y'_la `y`y'_la' text(`=`tick' + `y`y'_step'/100' `y`y'_x' " ", yaxis(`y') place(`y`y'_place') just(`y`y'_just') size(`y`y'_label_size') box lwidth(none) color(`y`y'_label_color') width(`y`y'_width') height(4) bcolor(white%70))
      *label
      local y`y'_la `y`y'_la' text(`tick' `y`y'_x' "{`y`y'_label_style':`y`y'_sign'`=string(`tick', "`y`y'_fmt'")'`y`y'_unit'}", yaxis(`y') place(`y`y'_place') just(`y`y'_just') size(`y`y'_label_size') box lwidth(none) color(`y`y'_label_color') bcolor(none))
    }
    else {
      local y`y'_la `y`y'_la' text(`tick' `y`y'_x' "{`y`y'_label_style':`y`y'_sign'`=string(`tick', "`y`y'_fmt'")'}", yaxis(`y') place(`y`y'_place') just(`y`y'_just') size(`y`y'_label_size') box lwidth(none) color(`y`y'_label_color') bcolor(none))
    }
  }
  
  * y-axis titles, horizontal orientation (build manually)
  local y 1
  local y1title text(`=`y`y'_tick_lo' + `y`y'_step'*(`y_ticks' - 1 + `y_margin_hi')' `xlo' "{`y`y'_title_style':Medical spending}"
  local y1title `y`y'title', yaxis(`y') place(`y`y'_place') just(`y`y'_just') box lwidth(none) color(`y`y'_title_color') bcolor(white) margin(b=+.5))
  local y1title `y`y'title' ytitle(" ", axis(`y') margin(l=-8))
  
  * y-axis settings
  local y_axis1 `y1_label' `y1_la' `y1_scale' `y1title'
  
  * Plot lines
  local y1 B
  local xvar x
  local line_y1_lfit "(lfit B x if UR==0, col(orange) lw(thick) range(200 800))"
  local line_y1_scat "(scatter B x if UR==0 & inrange(x, 201, 799),msym(S) mcolor(black) msize(*.8))"
  local line_y2_lfit "(lfit B x if UR==1, lpat(dash) lcol(orange) lw(thick) range(200 800))"
  local line_y2_scat "(scatter B x if UR==1 & inrange(x, 201, 799),msym(Oh) mcolor(black) msize(*1) mlwidth(medthick)) "
  local ci_area_0 (rarea Blb Bub `xvar' if UR==0, sort col(gs8%30) lw(none))
  local ci_area_1 (rarea Blb Bub `xvar' if UR==1, sort col(gs8%30) lw(none))
  
  * Label line_y1
  local text_1 text(12.51 500 "{`text_style':Mean unemployment}"
  local text_1 `text_1', yaxis(1) place(ne) just(left) box lwidth(none) width(45) color(`text_color') bcolor(white%90))
  local xref = `xvar'[4]
  sum `y1' in 4 
  local pci_1 (pci `=r(mean)' `=`xref'' `=r(mean) + .7' `=`xref' + 4.5', pstyle(p5))
  
  * Label line_y2
  local text_2 text(13 550 "{`text_style':+1 p.p. unemployment}"
  local text_2 `text_2', yaxis(1) place(ne) just(left) box lwidth(none) width(45) color(`text_color') bcolor(white%90))
  local xref = `xvar'[4]
  sum `y1' in 4 
  local pci_2 (pci `=r(mean)' `=`xref'' `=r(mean) + 0.07' `=`xref' + 0.45', pstyle(p5))
  
  grstyle set graphsize 3in 3in
  local legend legend(off)
  local region graphregion(margin(t=5 b=-1))
  local ws "placement(w) size(*1)"
  local es "placement(e) size(*1)"
 
  local labels          text(`=B[5]-0.5' `=x[5]-10' "{`text_style':age 49, mean}" "{`text_style':unemployment}", placement(s)) 
  local pci_labels      (pci `=B[5]' `=x[5]' `=B[5]-0.5' `=x[5]-10', pstyle(p5))

  local labels `labels' text(`=B[6]+0.5' `=x[6]+10' "{`text_style':age 49, +1 p.p.}" "{`text_style':unemployment}", placement(n))
  local pci_labels `pci_labels' (pci `=B[6]' `=x[6]' `=B[6]+0.5' `=x[6]+10', pstyle(p5))

  local labels `labels' text(`=B[7]-0.5' `=x[7]-10' "{`text_style':age 50, mean}" "{`text_style':unemployment}", placement(s)) 
  local pci_labels `pci_labels' (pci `=B[7]' `=x[7]' `=B[7]-0.5' `=x[7]-10', pstyle(p5))


  local labels `labels' text(`=B[8]+0.5' `=x[8]+10' "{`text_style':age 50, +1 p.p.}" "{`text_style':unemployment}", placement(n))
  local pci_labels `pci_labels' (pci `=B[8]' `=x[8]' `=B[8]+0.5' `=x[8]+10', pstyle(p5))

  twoway `ci_area_0' `ci_area_1' `line_y1_lfit' `line_y2_lfit' `line_y1_scat' `line_y2_scat' `pci_labels', `x_axis' `y_axis1' `legend' `region' `labels' graphregion(color(white))  bgcolor(white) 
  
  * Export graph
  cap mkdir "$SSDIMed/results/figures"
  graph export "$SSDIMed/results/figures/estimatedmodel_benefits4950_male`s'.pdf", as(pdf) fontface("Source Sans Pro") fontdir("$SSDIMed/data/raw/fonts") replace
}



